A3

(Zhu et al. 2008) (Xie 2023) (Robinson, McCarthy, and Smyth 2010) (Morgan 2022) (Durinck et al. 2009) (Gu et al. 2014) (Gu, Eils, and Schlesner 2016) (Wickham 2016) (Slowikowski 2023) (Raudvere et al. 2019)

if (!requireNamespace("GEOmetadb", quietly = TRUE))
  BiocManager::install("GEOmetadb")

if (!requireNamespace("knitr", quietly = TRUE))
  install.packages("knitr")

if (!require("edgeR", quietly = TRUE))
  BiocManager::install("edgeR")

if (!requireNamespace("BiocManager", quietly = TRUE))
  install.packages("BiocManager")

if (!requireNamespace("biomaRt", quietly = TRUE))
  BiocManager::install("biomaRt")

if (!requireNamespace("ComplexHeatmap", quietly = TRUE))
  BiocManager::install("ComplexHeatmap")

if (!requireNamespace("circlize", quietly = TRUE))
  BiocManager::install("circlize")

if (!requireNamespace("ggplot2", quietly = TRUE))
  install.packages("ggplot2")

if (!requireNamespace("ggrepel", quietly = TRUE))
  install.packages("ggrepel")

Calling required packages

library(BiocManager)
library(GEOmetadb)
library(knitr)
library(edgeR)
library(biomaRt)
library(ComplexHeatmap)
library(circlize)
library(ggplot2)
library(ggrepel)

Download the data

GSE152074 raw data supplementary file downloaded

sfiles = getGEOSuppFiles('GSE152075')

trying URL 'https://ftp.ncbi.nlm.nih.gov/geo/series/GSE152nnn/GSE152075/suppl//GSE152075_raw_counts_GEO.txt.gz?tool=geoquery'
Content type 'application/x-gzip' length 7121083 bytes (6.8 MB)
downloaded 6.8 MB

fnames = rownames(sfiles)
# there is only one supplemental file
readData = read.table(fnames[1],header=TRUE, check.names = TRUE)

Data

Data from (Lieberman et al. 2020).

kable(readData[1:5, 1:5], type = "html", row.names = TRUE)

	POS_001	POS_002	POS_003	POS_004	POS_005
A1BG	0	1	0	0	18
A1CF	0	0	2	0	0
A2M	69	36	84	42	83
A2ML1	2	0	0	0	3
A2MP1	0	0	0	0	21

Table 1: Original data contains HGNC annotation as row names. Column names have prefixes before their identifier number as either POS or NEG. Corresponding to either positive for COVID19 or negative.

Assess

Add 1 to all values of data so later on when conducting log2(cpm) we can avoid negative infinity values. (Advised by Professor Isserlin)

readData <- readData + 1

Setting first column as gene id for future format purposes

#Place rownames in first column for future format purposes
inter <- data.frame("HUGO" = rownames(readData))
geneData <- cbind(inter$HUGO, readData)
colnames(geneData)[1] <- "HUGO"

Clean

Remove any outliers that does not have at least 2 read per million in n of the samples. We set this as 2 since we add 1 to all of our dataset in the beginning of the code to have better plots. Denoting n as the smallest group of replicates which is the control group of 53. Using n = 53 conduct the removal of low counts.

#translate out counts into counts per millison using 
#the edgeR package function cpm
cpms = cpm(geneData[,2:485])
rownames(cpms) <- geneData[,1]
# get rid of low counts
keep = rowSums(cpms >2) >=53
geneData_exp_filtered = geneData[keep,]

Remove version numbers if they exists on gene id(HUGO) column. This makes it easier for mapping later on.

geneData_exp_filtered[,1] <- gsub("\\.[0-9]", "", geneData_exp_filtered[,1])

Map

#Mapping the name using biomatr
# list available gene annotation databases
bio <- useMart("ensembl", dataset = "hsapiens_gene_ensembl")
conversion_stash <- "geneMapping.rds"
if(file.exists(conversion_stash)){
  geneMapping <- readRDS(conversion_stash)
} else{
# convert column of gene IDs to Hugo symbols
geneMapping <- getBM(attributes = c("ensembl_gene_id", "hgnc_symbol"),
                     mart = bio,
                     filters = "hgnc_symbol",
                     values = geneData_exp_filtered[,1])
saveRDS(geneMapping, conversion_stash)
}

Combine the mapped gene data to original data

#Merge the data
mergedData <- merge(geneData_exp_filtered, geneMapping, by.x = 1, by.y = 2)
#remove duplicate rows in the gene data
mergedDataNoDup <- mergedData[!duplicated(mergedData[,1:485]),]

Apply Normalization

Randomly sample data to reduce the size of sample. Original sample is too large leading to computation errors due to the limitation of author’s computer.

set.seed(12345)
randomSamplePOS <- sample(mergedDataNoDup[2:431], 25)
randomSampleNEG <- sample(mergedDataNoDup[432:485], 25)
randomSample <- cbind(randomSamplePOS,randomSampleNEG, mergedDataNoDup$ensembl_gene_id, mergedDataNoDup$HUGO)

Define groups to use in normalization

samples <- data.frame(lapply(colnames(randomSample[1:50]), 
        FUN=function(x){unlist(strsplit(x, 
                        split = "_"))[c(2,1)]}))
colnames(samples) <- colnames(randomSample[1:50])
rownames(samples) <- c("patients","cell_type")
samples <- data.frame(t(samples))

Applying TMM to data

filtered_data_matrix <- as.matrix(randomSample[1:50])
rownames(filtered_data_matrix) <- randomSample$`mergedDataNoDup$ensembl_gene_id`
d = DGEList(counts=filtered_data_matrix, group=samples$cell_type)

d = calcNormFactors(d)

normalized_counts <- cpm(d)
#add columns of ensembl and hgnc id

normalized_count_data = data.frame(normalized_counts)
normalized_count_data$ensembl_gene_id <- mergedDataNoDup$ensembl_gene_id
normalized_count_data$hgnc_symbol <- mergedDataNoDup$HUGO


#This is a duplicate ensembl id that is giving errors when running code.
normalized_count_data <- normalized_count_data[-c(1902),]

model_design <- model.matrix(~samples$cell_type+0)
d <- estimateDisp(d, model_design)

Differential Gene Expression

LIMMA

p-value calculation using LIMMA

model_design <- model.matrix(~ samples$cell_type )

expressionMatrix <- as.matrix(normalized_count_data[,1:50])
rownames(expressionMatrix) <- 
  normalized_count_data$ensembl_gene_id
colnames(expressionMatrix) <- 
  colnames(normalized_count_data)[1:50]
minimalSet <- ExpressionSet(assayData=expressionMatrix)

Taking into account Patient variability

model_design_pat <- model.matrix(
  ~ samples$patients + samples$cell_type)

fit_pat <- lmFit(minimalSet, model_design_pat)

fit2_pat <- eBayes(fit_pat,trend=TRUE)

topfit_pat <- topTable(fit2_pat, 
                   coef=ncol(model_design_pat),
                   adjust.method = "BH",
                   number = nrow(expressionMatrix))
#merge hgnc names to topfit table
output_hits_pat <- merge(normalized_count_data[,51:52],
                         topfit_pat,by.y=0,by.x=1,all.y=TRUE)
#sort by pvalue
output_hits_pat <- output_hits_pat[order(output_hits_pat$P.Value),]

length(which(output_hits_pat$P.Value < 0.05))

[1] 1325

length(which(output_hits_pat$adj.P.Val < 0.05))

[1] 0

QLF

d = DGEList(counts=filtered_data_matrix, group=samples$cell_type)

d <- estimateDisp(d, model_design_pat)

fit <- glmQLFit(d, model_design_pat)

qlf.pos_vs_neg <- glmQLFTest(fit, coef='samples$cell_typePOS')
kable(topTags(qlf.pos_vs_neg), type="html",row.names = FALSE)

logFC	logCPM	F	PValue	FDR
8.508844	3.630861	56.94465	0	0.00e+00
9.217790	3.372834	44.85154	0	2.00e-07
8.238509	3.548781	39.51653	0	1.60e-06
8.480522	3.438800	38.21139	0	2.40e-06
7.357389	3.446313	35.77654	0	6.60e-06
14.042290	8.460151	34.52607	0	9.50e-06
8.560409	2.932483	34.40267	0	9.50e-06
10.527753	7.518887	32.50111	0	2.20e-05
7.998298	3.187208	31.11275	0	3.99e-05
7.616290	3.501343	29.92408	0	6.62e-05

x
BH

x
samples$cell_typePOS

x
glm

P-values were corrected using Quasilikelihood method. Quasilikelihood is better becacuse it is tailored towards RNAseq data

qlf_output_hits <- topTags(qlf.pos_vs_neg,sort.by = "PValue",
                           n = nrow(normalized_count_data))
length(which(qlf_output_hits$table$PValue < 0.05))

[1] 1569

length(which(qlf_output_hits$table$FDR < 0.05))

[1] 276

Thresholded over-representation analysis

Write to file upregulated, and downregulated genes

Which ones are upregulated and downregulated

length(which(qlf_output_hits$table$PValue < 0.05 
             & qlf_output_hits$table$logFC > 0))

[1] 1424

length(which(qlf_output_hits$table$PValue < 0.05 
             & qlf_output_hits$table$logFC < 0))

[1] 145

qlf_output_hits_withgn <- merge(randomSample[,51:52],qlf_output_hits, by.x=1, by.y = 0)
#number higher the lower the pvalue, and if it is upregulated number is positive, and negative for downregulated
qlf_output_hits_withgn[,"rank"] <- -log(qlf_output_hits_withgn$PValue,base =10) * sign(qlf_output_hits_withgn$logFC)
qlf_output_hits_withgn <- qlf_output_hits_withgn[order(qlf_output_hits_withgn$rank),]
upregulated_genes <- qlf_output_hits_withgn$`mergedDataNoDup$HUGO`[
  which(qlf_output_hits_withgn$PValue < 0.05 
             & qlf_output_hits_withgn$logFC > 0)]
downregulated_genes <- qlf_output_hits_withgn$`mergedDataNoDup$HUGO`[
  which(qlf_output_hits_withgn$PValue < 0.05 
             & qlf_output_hits_withgn$logFC < 0)]
write.table(x=upregulated_genes,
            file=file.path("data","upregulated_genes.txt"),sep = "\t",
            row.names = FALSE,col.names = FALSE,quote = FALSE)
write.table(x=downregulated_genes,
            file=file.path("data","downregulated_genes.txt"),sep = "\t",
            row.names = FALSE,col.names = FALSE,quote = FALSE)
write.table(x=data.frame(genename= qlf_output_hits_withgn$`mergedDataNoDup$HUGO`,F_stat= qlf_output_hits_withgn$rank),
            file=file.path("data","ranked_genelist.txt"),sep = "\t",
            row.names = FALSE,col.names = FALSE,quote = FALSE)

Non-thresholded Gene set Enrichment Analysis

gmt_url = "http://download.baderlab.org/EM_Genesets/current_release/Human/symbol/"
# list all the files on the server
filenames = RCurl::getURL(gmt_url)
tc = textConnection(filenames)
contents = readLines(tc)
close(tc)
# get the gmt that has all the pathways and does not include terms inferred
# from electronic annotations(IEA) start with gmt file that has pathways only
rx = gregexpr("(?<=<a href=\")(.*.GOBP_AllPathways_no_GO_iea.*.)(.gmt)(?=\">)", contents,
    perl = TRUE)
gmt_file = unlist(regmatches(contents, rx))
download.file(paste(gmt_url, gmt_file, sep = ""), "bader_lab.gmt")

trying URL 'http://download.baderlab.org/EM_Genesets/current_release/Human/symbol/Human_GOBP_AllPathways_no_GO_iea_April_02_2023_symbol.gmt'
Content type 'unknown' length 8507005 bytes (8.1 MB)
downloaded 8.1 MB

Conduct non-thresholded gene set enrichment analysis using the ranked set of genes from Assignment #2.  

1. What method did you use? What genesets did you use? Make sure to specify versions and cite your methods.  

I used the GSEA desktop application to run GSEA (Mootha et al. 2003).
I used the genesets from the bader lab extracted from the code above I used. Or it can be retrieved from here http://download.baderlab.org/EM_Genesets/current_release/Human/symbol/.

2. Summarize your enrichment results. 

SARs-CoV-2 positive samples:
4965 / 6074 gene sets are upregulated in phenotype na_pos
442 gene sets are significant at FDR < 25%
410 gene sets are significantly enriched at nominal pvalue < 1%
662 gene sets are significantly enriched at nominal pvalue < 5%
Top gene-set: HALLMARK_INTERFERON_ALPHA_RESPONSE%MSIGDBHALLMARK%HALLMARK_INTERFERON_ALPHA_RESPONSE
Number of genes in leading edge: 80
Top gene associated: CXCL11

SARs-CoV-2 negative samples:
1109 / 6074 gene sets are upregulated in phenotype na_neg
274 gene sets are significantly enriched at FDR < 25%
196 gene sets are significantly enriched at nominal pvalue < 1%
276 gene sets are significantly enriched at nominal pvalue < 5%
Top gene-set: ENERGY DERIVATION BY OXIDATION OF ORGANIC COMPOUNDS%GOBP%GO:0015980
Number of genes in leading edge: 165
Top gene associated: TEFM

3. How do these results compare to the results from the thresholded analysis in Assignment #2. Compare qualitatively. Is this a straight forward comparison? Why or why not?  

For the upregulated genes the top results are negative regulation of viral genome replication, negative regulation of viral process, response to type II interferon. downregulated we got cytoplasmic translation, positive regulation of respiratory burst, and intermediate filament-based process. For upregulated in thresholded and non-thresholded they align a bit together by having interferon related pathway results. Other then that they don’t seem to be similar. A common pathway in downregulated and and all of the genes for thresholded was a pathway related to cytoplasmic. However in non thresholded there were cytoplasmic related pathways but were very low in the list. It is not a straight forward comparison they both have different values as analysis. The thresholded could be more sensitive while the non-thresholded be more general.

Visualize your gene set Enrichment Analysis in Cytoscape  

1.  

Create an enrichment map - how many nodes and how many edges in the resulting map?  

379 Nodes
1511 Edges

What thresholds were used to create this map?

FDR q-value cutoff: 0.1
Analysis Type: GSEA
Node cutoff: 0.1
Edge cutoff: 0.375

Make sure to record all thresholds. Include a screenshot of your network prior to manual layout.  

Figure 0: Screenshot of network prior to manual layout

2. Annotate your network - what parameters did you use to annotate the network. If you are using the default parameters make sure to list them as well.

Used the auto annotate application in cytoscape.
Cluster algorithm: MCL Cluster
Label Column: GS_DESCR
Label Algorithm: WordCloud: Adjacent Words(default)
Max words per label: 3
Minimum word occurrence: 1
Adjacent word bonus: 8
Border Width: 3
Opacity: 20%
Font Scale: 20%

Figure 1: Annotated network

3. Make a publication ready figure - include this figure with proper legends in your notebook.

Figure 2: Collapsed

4. Collapse your network to a theme network. What are the major themes present in this analysis? Do they fit with the model? Are there any novel pathways or themes?  

Major themes present:

Up-regulated:

• GPCRS Rhodopsin ligand
  
• Chemokine Migration chemotaxis
  
• Cellular response necrosis
  
• type II interferon
  
• Migration monocyte chemotaxis
  
• Surface receptor pattern

Down-regulated:

• gluconeogenesis glycolysis
  
• superpathway warburg effect
  
• srp protein synthesis

Do they fit with the model?

There are major themes that coincide with the model present in the original paper. Such as the type II interferon, or surface receptor pattern is something that is common in COVID-19.

Novel pathways?
There exists some novel pathways in the network such as the foxo-mediated transcription of cell cycle gene. Pathways such as these there are not a lot of studies on them. Not just this one but out of the many pathways there are some novel interesting pathways to look at.

Interpretation and detailed view of results

The most important aspect of the analysis is relating your results back to the initial data and question.

1. Do the enrichment results support conclusions or mechanism discussed in the original paper? How do these results differ from the results you got from Assignment #2 thresholded methods.

This is a quotation taken from the original paper “SARS-CoV-2 induced a strong antiviral response with up-regulation of antiviral factors such as OAS1-3 and IFIT1-3 and T helper type 1 (Th1) chemokines CXCL9/10/11, as well as a reduction in transcription of ribosomal proteins” (Lieberman et al. 2020). They have found up-regulation in number of genes. The largest collapsed theme with 729 genes the GPCRS Rhodopsin ligand has two very large pathways named the GPCR ligand binding, and the CLASS A 1 (RHODOPSIN-LIKE RECEPTORS). Both of these gene-sets have at the top of their leading edge the CXCL 10/11/13 similar to the paper where they had up-regulation in the CXCL 9/10/11. In this sense the enrichment results support conclusions made in the original paper. Another example is the chemokine migration chemotaxis. In this major theme there are gene-sets involved in pathways such as the response to type II interferon, cellular response to type II interferon, . In the original paper they mention that as viral load increased the expression of interferon-responsive genes went up. Both of the two themes mentioned are up-regulated gene-sets. It seems that the enrichment results support the conclusion in the original paper.
Comparing from A2 threshold methods in the upregulated analysis using g:profiler we also had results such as response to type II interferon, and cellular response to type II interferon. This aligns with our GSEA results. The significant pathways seem to be aligning but the minor ones that have weak significance is where the major differences are betweeen t he Assignment 2 thresholded methods and Assignment 3.

For down regulated gene-sets for our network results from GSEA shows a major theme called srp protein synthesis. In this we have genesets such as the cytoplasmic translation. This gene-set is the largest gene-set inside srp protein synthesis and it is also the top result for down regulated genes gene-set for g:profiler using GO annotation.

Comparing the model in original paper, Assignment 2 results, and Assignment 3 results the significant pathways coincide but the less significant pathways is where the change occurs.

2. Can you find evidence, i.e. publications, to support some of the results that you see. How does this evidence support your result?

In a paper they state that SARS-CoV-2 infected thyroid gland activated the interferon pathways aligning with our results in the up-regulation (Poma et al. 2021). Another paper explains the down regulation of cytoplasmic translation. SARS-CoV-2 suppresses host protein translation leading to a down regulation of cytoplasmic translation pathways (Zhang et al. 2022). In both papers we can see supporting evidence for the results we got from the GSEA in upregulated and downregulated.

1. Add a post analysis to your main network using specific transcription factors, microRNAs or drugs. Include the reason why you chose the specific miRs, TFs or drugs (i.e publications indicating that they might be related to your model). What does this post analysis show?

Added a geneset for drugs from the bader lab. From this link (http://download.baderlab.org/EM_Genesets/current_release/Human/symbol/DrugTargets/Human_DrugBank_all_symbol.gmt).

I chose drugs because COVID-19 as of 2023 has no exact cure for the disease since the virus is constantly evolving. However having drug targets, and drugs to treat COVID-19 could be a very effective way of defeating the virus if we can reduce the risks of severity in COVID-19.

Conducted a exploratory analysis, and known signature analysis.
The exploratory analysis showed me the significant genesets and these are the ones I chose to focus on.
Mann-Whitney(Two-sided): 0.05

The most significant drug targets were Artenimol, Anisomycin, (S)-3-Phenyllactic Acid, Puromycin, NADH.

Using the Overlap is X percent of Sig gs: largest overlap drug targets is NADH, Aripiprazole, Loxapine, Ziprasidone, Aripiprazole Lauroxil

Hypergeometric Test: 0.25 The best test values are NADH, Aripiprazole, Loxapine, Ziprasidone, CYT997

In all 3 analysis we see that NADH is existent. NADH is the reduced form of NAD+. In the drugbank system they say that NADH has been useful in treating Parikinson’s disease. NADH is a significant nutrient in other diseases but how about for COVID-19. Multiple papers shows that NAD+ deficiency might be one of the main causes of the disease severity in COVID-19 (Miller, Wentzel, and Richards 2020). This aligns with the significant gene sets that are seen in the cytoscape analysis.

The figure is too large so none of the nodes are actually visible in cytoscape unless zoomed in.  Since it is difficult to view in network view. Here are the table views and just a few drug targets zoomed in for the analysis done.  

Figure 5: Drug Targets zoomed in Mann-Whitney network

Figure 6: Drug Targets Mann-Whitney table view

Figure 9: Drug Targets zoomed in overlap

Figure 10: Drug Targets table overlap

Compilation

This code compiles with docker as of 04/04/2023.

References

Durinck, Steffen, Paul T. Spellman, Ewan Birney, and Wolfgang Huber. 2009. “Mapping Identifiers for the Integration of Genomic Datasets with the r/Bioconductor Package biomaRt.” Nature Protocols 4: 1184–91.

Gu, Zuguang, Roland Eils, and Matthias Schlesner. 2016. “Complex Heatmaps Reveal Patterns and Correlations in Multidimensional Genomic Data.” Bioinformatics. https://doi.org/10.1093/bioinformatics/btw313.

Gu, Zuguang, Lei Gu, Roland Eils, Matthias Schlesner, and Benedikt Brors. 2014. “Circlize Implements and Enhances Circular Visualization in r.” Bioinformatics 30: 2811–12.

Lieberman, Nicole AP, Vikas Peddu, Hong Xie, Lasata Shrestha, Meei-Li Huang, Megan C Mears, Maria N Cajimat, et al. 2020. “In Vivo Antiviral Host Transcriptional Response to SARS-CoV-2 by Viral Load, Sex, and Age.” PLoS Biology 18 (9): e3000849.

Miller, R., A. R. Wentzel, and G. A. Richards. 2020. “COVID-19: NAD+ Deficiency May Predispose the Aged, Obese and Type2 Diabetics to Mortality Through Its Effect on SIRT1 Activity.” Medical Hypotheses 144: 110044. https://doi.org/https://doi.org/10.1016/j.mehy.2020.110044.

Mootha, VK, CM Lindgren, KF Eriksson, A Subramanian, S Sihag, J Lehar, P Puigserver, et al. 2003. “PGC-1α-Responsive Genes Involved in Oxidative Phosphorylation Are Coordinately Downregulated in Human Diabetes.” Nature Genetics 34 (3): 267–73. https://doi.org/10.1038/ng1180.

Morgan, Martin. 2022. BiocManager: Access the Bioconductor Project Package Repository. https://CRAN.R-project.org/package=BiocManager.

Poma, Anna Maria, Alessandro Basolo, Daniele Bonuccelli, Agnese Proietti, Elena Macerola, Clara Ugolini, Ludovica Torregrossa, et al. 2021. “Activation of Type i and Type II Interferon Signaling in SARS-CoV-2-Positive Thyroid Tissue of Patients Dying from COVID-19.” Thyroid : Official Journal of the American Thyroid Association 31 (12): 1766–75. https://doi.org/10.1089/thy.2021.0345.

Raudvere, Uku, Liis Kolberg, Ivan Kuzmin, Tambet Arak, Priit Adler, Hedi Peterson, and Jaak Vilo. 2019. “G:profiler: A Web Server for Functional Enrichment Analysis and Conversions of Gene Lists.” Nucleic Acids Research. https://doi.org/10.1093/nar/gkz369.

Robinson, Mark D, Davis J McCarthy, and Gordon K Smyth. 2010. “edgeR: A Bioconductor Package for Differential Expression Analysis of Digital Gene Expression Data.” Bioinformatics 26 (1): 139–40. https://doi.org/10.1093/bioinformatics/btp616.

Slowikowski, Kamil. 2023. Ggrepel: Automatically Position Non-Overlapping Text Labels with ’Ggplot2’. https://CRAN.R-project.org/package=ggrepel.

Wickham, Hadley. 2016. Ggplot2: Elegant Graphics for Data Analysis. Springer-Verlag New York. https://ggplot2.tidyverse.org.

Xie, Yihui. 2023. Knitr: A General-Purpose Package for Dynamic Report Generation in r. https://yihui.org/knitr/.

Zhang, Dongdong, Longlong Zhu, Yu Wang, Pengfei Li, and Yu Gao. 2022. “Translational Control of COVID-19 and Its Therapeutic Implication.” Frontiers in Immunology 13: 857490. https://doi.org/10.3389/fimmu.2022.857490.

Zhu, Yuelin, Sean Davis, Robert Stephens, Paul S. Meltzer, and Yidong Chen. 2008. “GEOmetadb: Powerful Alternative Search Engine for the Gene Expression Omnibus.” Bioinformatics (Oxford, England) 24 (23): 2798–2800. https://doi.org/10.1093/bioinformatics/btn520.